function grav_run_stage4(fname3)
  % make hydro files callable from the current directory
  thisdir = fileparts(mfilename('fullpath'))
  addpath(fullfile(thisdir, '../hydro/'))


  datestring = datestr(clock, 'yy-mm-dd-HHMM');


  Nz = 32;

  load(fname3, 'jobname', 'params_hi', 'mep_tsteps', 'mep_time', 'mep_dt');

  Nx = params_hi.Nx;

  params_grav = grav_params(params_hi, Nz);
  params_grav.metricLinearFnInterpolation = true;
  params_grav.plotMetric = true;

  grav_init(params_grav);

  %*********************************************
  %  load hydro solution 
  %*********************************************
  load(fname3, 'uOut3', 'udtOut3', 'udtdtOut3');



  % vector of metric orders to evaluate
  %metricorder = [0 1 2];
  metricorder = [1];


  % read off the number of evaluation points
  Nmep = length(mep_time);

  % read off the differencing order
  tord = size(mep_tsteps, 2) - 1;

  for i=1:Nmep
    params_grav.t  = mep_time(i);
    params_grav.dt = mep_dt(i);

    params_grav

    etrh   = zeros(length(metricorder), params_grav.Nx, params_grav.Ny);
    detrb  = zeros(length(metricorder), params_grav.Nx, params_grav.Ny);
    etri   = zeros(length(metricorder), params_grav.Nx, params_grav.Ny);
    etr2i  = zeros(length(metricorder), params_grav.Nx, params_grav.Ny);

    % prepare hydro solution
    u     = squeeze(uOut3(i, :, :));
    udt   = squeeze(udtOut3(i, :, :));
    udtdt = squeeze(udtdtOut3(i, :, :));

    for j=1:length(metricorder)
      params_grav.metricOrder = metricorder(j);
      [g gt gtt] = grav_metric(params_grav, u, udt, udtdt);
      fname = sprintf('%s-%s-s4-o%d-t%d', datestring, jobname, metricorder(j), round(mep_time(i)));
      
      save(fname, 'jobname', 'params_grav', 'g', 'gt', 'gtt', '-v7.3');
    end
  end

end
%       tic
%       [etrh(j,:,:,:), detrb(j,:,:,:), etri(j,:,:,:), etr2i(j,:,:,:)] = grav_grtensors(params_grav, g, gt, gtt);
% 
%       disp(sprintf('evaluated gr tensors in %.2f sec', toc))
%     end
%     fname = sprintf('%s-%s-gr-t%d', datestring, jobname, round(mep_time(i)));
%     save(fname, 'jobname', 'params_grav', 'etrh', 'detrb', 'etri', 'etr2i', '-v7.3');
%   end
% 
%   
% 
% end


















